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The Difference PDF of 21-cm Fluctuations: A Powerful 
Statistical Tool for Probing Cosmic Reionization 
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ABSTRACT 

A new generation of radio telescopes are currently being built with the goal of tracing 
the cosmic distribution of atomic hydrogen at redshifts 6-15 through its 21-cm line. The 
observations will probe the large-scale brightness fluctuations sourced by ionization 
fluctuations during cosmic reionization. Since detailed maps will be difficult to extract 
due to noise and foreground emission, efforts have focused on a statistical detection of 
the 21-cm fluctuations. During cosmic reionization, these fluctuations are highly non- 
Gaussian and thus more information can be extracted than just the one-dimensional 
function that is usually considered, i.e., the correlation function. We calculate a two- 
dimensional function that if measured observationally would allow a more thorough 
investigation of the properties of the underlying ionizing sources. This function is 
the probability distribution function (PDF) of the difference in the 21-cm brightness 
temperature between two points, as a function of the separation between the points. 
While the standard correlation function is determined by a complicated mixture of 
contributions from density and ionization fluctuations, we show that the difference 
PDF holds the key to separately measuring the statistical properties of the ionized 
regions. 
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1 INTRODUCTION 

The earliest generations of stars are thought to have trans- 
formed the universe from darkness to light and to have 
reionized and heated the intergalactic medium. Knowing 
how the reionization process happened is a primary goal 
of cosmologists, because this would tell us when the early 
stars formed and in what kinds of galaxies. The cluster- 
ing of these galaxies is particularly interesting since it is 
driven by large-scale de nsity fluctuations in the dark matter 
<|Barkana fe Loebl l2004h . While the distribution of neutral 
hydrogen during reionization can in principle be measured 
from maps of 21-cm emission by neutral hydrogen, upcoming 
experiments such as the Mileura Widefield ArrajQ and the 
Low Frequency ArrajQ are expected to be able to detect ion- 
ization fluctuations only statistically (for rev iews see, e.g., 
iFurlanetto et alJl200d : lBarkana fc Loebll2007ft . 

Studies of statistics of the 21-cm fluctuations have 
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focused on the two-point correlation function (or power 
spectrum) of the 21-cm brightness temperature. This is 
true both for analytical and numerical studies and anal- 
yses of the expected sensitivity o f the new experiment s 
jBowman. Morales, fc Hewittl l200rj ; iMcQuinn et all 120061 ). 
The power spectrum is the natural statistic at very high red- 
shifts, as it contains all the available statistical information 
as long as Gaussian primordial density fluctuations drive the 
21-cm fluctuations. However, during reionization the hydro- 
gen distribution is a highly non-linear function of the distri- 
bution of the underlying ionizing sources. This follows most 
simply from the fact that the H I fraction is constrained to 
vary between and 1, and this range is fully covered in any 
scenario driven by stars, in which the intergalactic medium 
is sharply divided betwe en H I and H II region s. The re- 
sulting non-Gaussianity (Bh aradwai fc A"lHl2005l ) raises the 
possibility of using complementary statistics to measuring 
additional informati on that is not directly derivable from 
the power spectrum l|Saivad-Ali et al1 [2006). 

Numerical simulations have recently begun to reach the 
large scales (of order 100 Mpc) neede d to capture the evo- 
lution of the IGM during reionization (Mcll ema et aljfeood ; 
IZahn et al Jl2007h . These simulations account accurately for 
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gravitational evolution on a wide range of scales but still 
crudely for gas dynamics, star formation, and the radiative 
transf er of ionizing photons. Analytically, iFurlanetto et al] 
used the statistics of a random walk with a linear bar- 
rier to model the H II bubble size distribution during the 
reionization epoch. Schematic approxima tions were devel- 
oped for the two-point corr elation functi on (IFurlanetto et all 
l2004lMcQuinn et alj2005( ). but recentlv lBarkanal (|2007f t de- 
veloped an accurate, self-consistent analyti cal expression for 
the fu ll two-point distribution within the IFurlanetto et al] 
model, and in particular used it to calculate the 21- 
cm correlation function. 

Noting the expected non -Gaussianity and the i mpor- 
tance of additional statistics, IFurlanetto et al] l|2004 ) also 
calculated the one-point probability distribution function 
(PDF) of the 21-cm brightness temperature Tt at a point. 
The PDF has begun to be explored in numerical simula - 
tions as well i|Ciardi fc Madaull2003l ; iMellema et al.ll2006ft . 
Some of the additional informati on available in the PDF 
can be capture d by the skewness (IWvithe fc Morales! 120071 ') 
or bispectrum (Sai vad-Ali et alll2006l ) statistics. Both the 
correlation function and the PDF are functions of a single 
variable (at each redshift): the two-point correlation func- 
tion is a function of separation, and the PDF is a function 
of T(,. It is possible to create a two-dimensional function by 
calculating the one-point PDF as a function of smoothing 
scale (or pixel size), but this quantity is difficult to inter- 
pret since it is not simply related to the 21-cm correlation 
function or to the ionization statistics. 

In this paper we consider a two-dimensional function 
that generalizes both the one-point PDF and the correla- 
tion function and yields additional information beyond those 
statistics. In particular, the variance of this new statistic is 
simply related to the 21-cm correlation function that is usu- 
ally considered. This function is the PDF of the difference 
AT;, = Tb,i — Tb,2 of the 21-cm brightness temperatures Tb 
at two points. We present in the next section our analytical 
model for predicting the difference PDF; its precise relation 
to the two-point correlation function is presented in sec- 
tion [231 We present illustrative predictions of the difference 
PDF in section 3, where in section 3.2 we emphasize that 
it can be used to separately measure ionization correlations. 
We summarize our conclusions in section 4. 



2 MODEL 

Analytical approaches to galaxy formation and reionization 
are based on the mathematical problem of random walks 
with barriers. The statistics of a single random walk can be 
used to calculate various one-point distributions; in partic- 
ular, the statistics of a random walk with a linear barrier 
can be used to calculate th e distribution of ionize d bub- 
ble sizes during reionization (jFurlanetto et all 120041 ). How- 
ever, to calculate the correlation function and other two- 
point distributions requires us to solve for the simultane- 
ous evolution of two correlated random wal ks at two dif- 
ferent points. IScannapieco fc Barkanal l|2002h found an ap- 
proximate but quite accurate analytical solution in the case 
of constant barriers and used it to calculate the joint, bi- 
variate mass function o f hal os forming at two redshifts; 
IScannapieco fc Thackerl (|2005T ) showed that this solution de- 



scribes well the two-point correlation function of halos in 
numerical simulations, particularly when expressed in La- 
grangian coordin ates (i.e.. in ter ms of the initial comoving 
halo separation). iBarkanal (|2007l ) generalized the two-point 
solution to the case of linear barriers, and applied it to cal- 
culate the correlation function of cosmological 21-cm fluc- 
tuations during reionization. In this section we first review 
the basic setup of the two-barrier problem in the context 
of reioni z ation . We then briefly summarize the solution of 
IBarkanal (|2007T l, except that we generalize it slightly to the 
case of measurements in the presence of additional smooth- 
ing (e.g., due to a limited instrumental resolution). We show 
how to apply this solution to calculate the difference PDF 
of 21-cm fluctuations. 

2.1 Reionization: basic setup 

The basic approa ch for using r andom walks with barriers in 
cosmology follows iBond et ail l|l99ll ). wh o used it to rederive 
and e xtend the halo formation model of iPress fc Schechterl 
(1 19741 ). In this approach we work with the linear overden- 
sity field <5(x, z) = p(x,z)/p(z) — 1, where x is a comoving 
position in space, z is the cosmological redshift and p is the 
mean value of the mass density p. In the linear regime, the 
overdensity grows in proportion to the linear growth fac- 
tor D{z) (defined relative to z = 0). This fact is used in 
order to extrapolate the linear density field to the present 
time, i.e., the initial density field at high redshift is extrapo- 
lated to the present by multiplication by the relative growth 
factor. We adopt this view, and throughout this paper quan- 
tities such as S and the power spectrum P(k) refer to their 
values linearly-extrapolated to the present. In each applica- 
tion there is in addition a barrier that signifies the critical 
value which the linearly-extrapolated S must reach in order 
to achieve some physical milestone on some scale. In this 
work the milestone corresponds to having a sufficient num- 
ber of galaxies within some region in order to fully reionize 
that same region. 

At a given z, we consider the smoothed density in a re- 
gion around a fixed point A in space. We begin by averaging 
over a large scale or, equivalently, by including only small co- 
moving wavenumbers k. We then average over smaller scales 
(i.e., include larger k) until we find the largest scale on which 
the averaged overdensity is higher than the barrier; in the 
application to reionization, we then assume that the point 
A belongs to an H II bubble of this size. Mathematically, if 
the initial density field is a Gaussian random field and the 
smoothing is done using sharp fc-space filters, then the value 
of the smoothed S undergoes a random walk as the cutoff 
value of k is increased. Instead of using k, we adopt the 
(linearly-extrapolated) variance S of density fluctuations as 
the independent variable. While the solutions are derived 
in reference to sharp fc-space smoothing, we follow the tra- 
ditional extended Press-Schechter approach and substitute 
real-space quantities in the final formulas. In particular, S 
is calculated as the variance of the appropriate mass M en- 
closed in a spatial sphere of comoving radius r. 

We apply mathematical random-walk statistics to the 
distribut ion of H II regions dur ing reionization using the 
model of IFurlanetto et ail (|2004f ). According to this model, 
a given point A is contained within a bubble of size given 
by the largest surrounding spherical region that contains 
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enough ionizing sources to fully reionize itself. If we ignore 
recombinations, then the ionized fraction in a region is given 
by C/coil, where / co ii is the collapse fraction (i.e., the gas 
fraction in galactic halos) and £ is the overall efficiency fac- 
tor, which is the number of ionizing photons that escape 
from galactic halos per hydrogen atom (or ion) contained in 
these halos. This simple version of the model remains ap- 
proximately valid even with recombinations if the number 
of recombinations per hydrogen atom in the IGM is roughly 
uniform; in this case, the resulting reduction of the ionized 
fraction by a constant factor can be incorporated into the 
value of (. 

In the extended Press-Schechter model l|Bond et al.l 

Il99ll ), in a region containing a mass corresponding to vari- 
ance S, 



/coii = erfc 



S c (z) - 8 



S) 



(1) 



where Smi n is the variance corresponding to the minimum 
mass of a halo that hosts a galaxy, 8 is the mean den- 
sity fluctuation in the given region, and S c (z) is the crit- 
ical density for halo collapse at z. While this describes 
fluctuations in / co n well, the cosmic mean collapse frac- 
tion (and thus the overall evolution of reionization with 
redshift) is better desc ribed by the halo mass function of 
Sheth & Tormcn ( 1999) (with the updated parameters sug- 
gested by ISheth fc Tormenl (120021 )'). We thus use the lat- 
ter mean mass function and adjust / co ii in different re- 
gions in proportion to th e extended Press-Schechter formula; 
iBarkana fc Loebl |2004) suggested this hybrid prescription 
and showed that it fits a broad range of simulation results. 
With these assumptions, the exact ionized fraction x 1 in a 
region is given by 



•^exact i z > 5mm, 8, S) — £ = 

JPS 



\ v /2(S* mi „ - S) t 



(2) 



where /st and /ps are the cosmic mean collapse fractions 
according to the Sheth- Tormen and Press-Schechter models, 
respectively. 

The resulting condition x % — 1 for having an ionized 
bubble of a given size, written as a cond i tion fo r 8 vs. S, is 
of the same form as in iFurlanetto et al] (120041 ). at a given 
redshift, and thus (as they showe d) yields a linear barri er 
to a good approximation (see also IFurlanetto et al.l |2006)). 
We write the effective linear barrier in a general notation: 



S = v + fiS , 
where 



(3) 



8c{z) — V 25'min erfc 



C/ST 



^ = 



■. erfc 



VC/st/ 



(4) 



with erfc - denoting the inverse function of erfc. For con- 
sistency, in what follows we use a modified formula for the 
ionized fraction, replacing equation @ with the expression 
that corresponds to the linear approximation of the barrier: 



x\z, Smm, 5, S) = C erfc 



8 c (z) 



VSSmin (l - 2S~) 



(•») 



This replacement ensures that the ionized fraction varies 
from to I as 8 goes from —00 up to the barrier. We also 
denote the neutral fraction x n — 1 — x l and, in particular, 



X (2, tS'min, 8, 5*) 



X (z, jSmin, 8, iS*) 



(6) 



The approximation of the linear barrier is quite accurate as 
long as the maximum S that we consider is much smaller 
than Smi n , which is the case in the applications of the model 
in this paper, where the maximum S is set by the resolution 
of the upcoming experiments (see the next subsection). 

Because of some approximations in this model, the total 
ionized fraction as given by the model [see equation (|15p in 
the next subsection] comes out slightly different from the di- 
rect result for the mean global ionized fraction, (x l ) = C/st 
in terms of the cosmic mean coll apse fraction (Note that the 
model of IFurlanetto et all (|2004l ) suffers from a similar diffi- 
culty). To deal with this, we adopt the direct values of (x 1 ) 
versus redshift, and adjust £ within the model to an effective 
value of £ at each redshift that gives a model value of (x l ) 
that equals the desired one. This typically only requires an 
adjustment by a few percent or less. 



2.2 The 21-cm one-point PDF 

Before considering two-point functions, we fi rst calculate the 
21-cm PDF around one point by following IFurlanetto et all 
(2004), except that we obtain the ionized fraction from 
cq. ((5]l for consistency with the barrier, and we also apply 
a non-linear correction to the density. We denote the PDF 
itself by Pr b , and the cumulative probability distribution 
(CPD) J P Tb dT b by C Tb - 

During cosmic reionization, we assume that there are 
sufficient radiation backgrounds of X-rays and of Lya pho- 
tons so that the cosmic gas has been heated to well above 
the cosmic microwave background temperature and the 21- 
cm level occupations have come into equilibrium with the 
gas temperature. In this case, the observed 21-cm brightness 
temperature relative to the CMB is independent of the spin 
temperatur e and, for our assum ed cosmological parameters, 
is given by l|Madau et al.ll 19971 ) 



T b = T 6 * ; T b = 25 



1+z 



8 



mK 



(7) 



with ^ = x n [l + 8 L {z)\, where x n is the neutral hydrogen 
fraction and 8 L {z) = D(z)8 is the linear overdensity at 2 (as 
opposed to 8 which denotes the density linearly-extrapolated 
to redshift 0). Under these conditions, the 21-cm fluctuations 
are thus determined by fluctuations in 

In the model, x n is determined by the halo abundance, 
which is in turn determined by the statistics of the linear 
density field, and thus x n naturally falls within the correct 
range of to 1. However, ^ also depends on the actual den- 
sity, and the linear density S L (z) can take on unphysical val- 
ues below -1. While a full non-linear model would be difficult 
to solve, within the context of the model where statistics are 
averaged over spherical regions, we can make a simple, ap- 
proximate correction in order to get reasonable values for 
the actual, nonlinear density 8 NL (z). iMo fc White] (Il996l ) 
developed such an approximate formula for <5 L as a func- 
tion of (5 NL , based on spherical collapse (of overdensities) or 
spherical expansion (of voids). In particular, they incorpo- 
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rated the asymptotic limits of <5 NL — + — 1 and <5 NL — *• oo as 
well as the correct behavior near <5 NL = 0. Their formula is 
valid for the Einstein-de Sitter universe or more generally at 
high redshift (when the dynamical effect of the cosmological 
constant is negligible) . We similarly develop and use here an 
accurate approximation for the inverse function, 



1 + 5 NL (<5 L ) 



1.12431 



+ 



0.395223 (Sr 



1.35 

L n -1.72256 



-3/2 



(8) 



where Sf — 1.68647 is the critical collapse overdensity in an 
Einstein-de Sitter cosmology. Thus, the expression we use 
for ^ is 



* = x n {l + S NL [D(z)S]} 



(9) 



If we denote by P* the one-point PDF of at a point, 
then this is related to the 21-cm one-point PDF defined 
above by 



P Tb {T b ) = ^Pv 

lb 



(10) 



Also, the corresponding cumulative probability distributions 
are equal. We also assume that the PDF is considered on 
a resolution scale r rcs (corresponding to a variance S Tes ), 
i.e., that the density and ionization states are averaged on 
the scale r las around the point being considered. In other 
words, we are really considering the density and ionization 
distributions in a region of size r res centered at a point. 

We must now consider the separate contributions to 
the PDF of ^ from two possible cases. First, if the point 
lies within a fully ionized region, then ^ is identically zero, 
so this case contributes a fo-function (Dirac delta func- 
tion) at 0, containing the total probability that the re- 
gion is ionized. This probability i s given by the quantity 
f>,iin(^, Mi ffres) in equation (15) in lBarkanal (|2007T l. as first 
derived by McQuin n et all (|2005l ). The second case is if the 
point lies within a region that is still partially neutral. In 
this case, VP is given by equation (JSJ, where for each value 
of 8 we use x n = x n (z, Smin, 5, S les ) from equation ||BJ. We 
then distribute the total probability of this case into various 
values of 9 using the conservation of probability, i.e., for 
each possible val ue of S, the prob ability Qu^is, \i, 5, S les ) dS 
[equation (14) i n Barkanal (|2007f ). again in agreement with 
iMcQuinn et"afl l|2005h ] contributes to P* ( ^) at the value 
of ^ that corresponds to 8. As noted by iFurlanetto et al] 
(|2004| ), there is a maximum value of $ , $ max , since \P — > 
in the two extreme limits, both when the density goes to 
zero (due to the density term in and when it goes up 
to the barrier (due to full ionization). Thus, two values of 
8 contribute to each value of ty, and P* (*&) is singular at 

*max- 



2.3 The 21-cm difference PDF 

We now consider two separate points, with the same as- 
sumptions as in the previous subsection. We wish to consider 
the PDF of the difference AT b = T M - T b>2 of the 21-cm 
brightness temperatures T b at two points (or, in fact, aver- 
aged over regions centered at each of the two points). We 
denote the PDF itself by P A T b , and the CPD / P A T b dAT b 



by Cat,,- If we denote by Pa* the PDF of the difference 
A$ = — ^2 1 between the values of $ in the two regions, 
then this is related to the 21-cm difference PDF by 



s"" (It) 



("I 



Also, the corresponding cumulative probability distributions 
are equal. 

We therefore wish to determine the PDF of the differ- 
ence as a function of the comoving distance d between 
two points being considered at redshift z. As before, we also 
assume that the PDF is considered on a resolution scale 
r rcs , i.e., we are really considering the joint density and ion- 
ization distributions of two regions of size r las ce ntered at 
two p oints separated by a distance d. The model of lBarkanal 
()2007h provides the probability that either or both of these 
regions lie completely within H II bubbles and, when the re- 
gions are not fully ionized, the model provides the correlated 
distributions of their average overdensities. 

We must now consider the separate contributions to the 
PDF of A>I> from three possible cases. First, if both points 
are within fully ionized regions, then A^ is identically zero, 
so this case contributes a <5_D-function at 0, containing the 
total probability that both regions are ionized. This proba- 
bility is given by th e quantity P> (v v, fi, fi, S ICS , S ICS , £res(d)) 
in equation (40) in IBarkanal (|2007) . where £ is the effective 
real-space cross-correlatio n between the d ensities of the two 
regions [see section 4.1 of lBarkanal (|2007h ]: 



kd 



(12) 



where W(x) is the Fourier transform of a spherical top-hat 
window function. The second case is where one of the re- 
gions is fully ionized, so, e.g., we assume that region 2 is in 
an H II bubble while region 1 is not, and then double the 
contribution in order to include the symmetric, opposite sit- 
uation. In this case, A* = *i = x {l + 8 Nlj [D{z)5i\) , 
where for each value of 8\ we use a;™ = x n (z, S m ln, Si, &es) 
from equation ((6]l . We then distribute the total probability of 
this case into various values of A^f using the conservation of 
probability, i.e., for each possible value of 8i, the probability 
/[fl|<>1 {v , y, MiM i Si, S les , S les , £iea(d)) dSi [equation (43) in 
IBarkanal (2007)] contributes to P A *(A*) dA* at the value 
of A^ that corresponds to Si. The final case is where both 
regions are not fully ionized. In this case, A^ = — 9 2 \ is 
a function of Si and 82, and the conservation of probabil- 
ity turns Q(i>, v, fi, 11 , 81,82 , S TCB , S TCS , £tes(d)) dSi dSi [equa- 
tion (36) in lBarkanal <|2007f )] into P A *(A*) dA* at the ap- 
propriate value of A<3/. 

The variance of the 21-cm PDF (for two points sepa- 
rated by a distance d) is 

(AT 6 2 > = 2[CAT b (0)-£AT b (d)] , (13) 

where the ordinary 21-cm correlation function is ^AT b (d) = 
(AT 6l iAT(,,2) - (AT M }(ATf,,2), and the f notation denotes 
averaging on the resolution scale r rGS . The variance can be 
calculated from the PDF using one of these expressions: 



(ATi) = / (AT b ) 2 PAT b {AT b )dAT b 



2 / A% [1 - C*AT b (AT b )} dAT b 



(14) 



The Difference PDF of 21-cm Fluctuations: A Powerful Statistical Tool for Probing Cosmic Reionization 5 



where we integrated by parts to get the second expression. 

To calculate the correlation function of Tb (or, equiv- 
alently, of without first calculating the PDF , we can 
calcul ate various expectation values using the IBarkanal 
(2007) solution, once again generalized to include a reso- 
lution/smoothing length. First, the mean ionized fraction in 
the model is 



(x % ) = F > A in (u, fl, SVcs) 

— I Ql\ n (v,^,8,S ICS )x{z,S m in,8,S T es)d8 



(15) 



where f>,ii n and Qu n are given in section 3 of IBarkanal 
(2007), while the mean $ at a point is 

{9) = l-{x i )+ (16) 

Ql in (l/, fl, 5, S T m) X n (z, S m in, 6, S rcs ) S NL (D(z)S) dS . 
S= — oo 

For two points separated by a distance d, 



d5i 



dS 2 



(17) 



Q{v, V., (JL, fl, 5l,52, Srcs, S x es, ^roa(rf)) 

x [l+8 Nh (D(z)8 1 )] [l + 8 NL (D(z)S 2 )] 

x x n (z, S min , ^1 , Srcs ) X (z, aSmm, $2 , Sres) , 

which is a generalization of equation (49) of Bar kanal (|2007l ) 
and reduces to that equation in the limit S las — > Smin- Also, 
in the limit d — > equation (|17p simplifies to 

(* 2 ) = / d8Q lin (u,^,S,Sr CS ) (18) 



x {[l + 8 Nh (D(z)8)] x n (z,S min ,S,S vcs )Y 



3 RESULTS 

3.1 The full 21-cm difference PDF 

In this subsection we use our model to predict the 21-cm 
difference PDF, plot it in full as a two-dimensional func- 
tion of the separation d and the brightness temperature dif- 
ference ATb of the two points, and explore its dependence 
on a number of the input parameters (including the red- 
shift and the resolution scale). In the following subsection 
we then show that important information can be extracted 
using gross features of the PDF that are insensitive to its 
detailed shape; in particular, this information can be used to 
cleanly separate out and measure statistics of the ionization 
field that otherwise would be mixed in and convolved with 
the density field within the usually-considered two-point cor- 
relation function. Throughout this section, we illustrate our 
predictions in a ACDM universe that includes dark matter, 
baryons, radiation, and a cosmological constant. We assume 
cosmological parameters that match the three year WMAP 
data together with weak lensing observations (jSpergel et al.l 
I2007D . namely fi m = 0.299, tt A = 0.701, fi 6 = 0.0478, 
h = 0.687, n = 0.953 and cr 8 = 0.826. 

Figure [T] shows an example of the 21-cm one-point PDF 
and CPD and the two-point 21-cm difference PDF and CPD. 
For the one-point function, there are two separate contri- 
butions; the case of a partly neutral region is shown as a 
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Figure 1. The one-point probability distribution function (PDF) 
and cumulative probability distribution (CPD) of the 21-cm 
brightness temperature (panels on the left), and the two-point 
difference PDF and CPD (panels on the right). We consider red- 
shift z = 8.26, when the mean global ionized fraction is (x l ) = 0.5. 
In our calculations we assume a constant ionization efficiency in 
all halos in which atomic cooling occurs (i.e., in all halos above 
a circular velocity V c = 16.5 km/s), set so that reionization ends 
at z = 6.5. The assumed resolution (i.e., the diameter of each 
region) is 1' in all panels. Note that in each case, the PDF has an 
additional <5d -function (not shown) that brings the value of the 
corresponding CPD to unity at AT;, = (or = 0). Also, the 
PDF becomes singular at the maximum value of ATj (or T;,). For 
the two-point distributions, the full PDF or CPD (solid curve) is 
the sum of contributions from both regions being neutral (short- 
dashed curve) and from one of them being neutral (long-dashed 
curve); the comoving separation between the two centers is set 
equal to 10 Mpc. 



function of Tb, while the case where the region is ionized 
contributes an additional <5.o-function to the PDF, or equiv- 
alently a step function to the CPD. In the CPD the size of 
this step function can be easily read off as the additional 
value needed to bring it up to unity at 21, = 0. Similarly, 
for the difference PDF and CPD there are three separate 
contributions; two of them - the cases of both regions being 
partly neutral or just one of them - are shown as functions 
of ATb, while the third case - with both regions fully ion- 
ized - contributes an additional OD-function to the PDF. 
Again, in the CPD the size of the step function equals the 
additional value needed to bring the CPD up to unity at 
ATb = 0. Another adv antage of the one - point or two-point 
CPDs, pointed out by iFurlanetto et al.l (|2004T > in the one- 
point case, is that the PDF becomes singular (in the model) 
at the maximum value of Tb while the CPD does not diverge. 
Thus, for these two reasons, we henceforth prefer to plot the 
CPD instead of the PDF. 

The model predicts characteristic shapes for the PDF 
and CPD during reionization. In particular, the one-point 
function cuts off at some maximum value of ATb (which 



6 R. Barkana and A. Loeb 



1 .2 r — i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 




<5Kz) 



Figure 2. The value of Vt 1 and x n in a region versus the mean 
linear 8 in that region at rcdshift z. We show the value of \E> as 
given by eq. J9} [solid curves] and of the mean neutral fraction 
x n in the region as given by eq. [dashed curves]. We consider 
redshifts when the mean global ionized fraction is (x 1 ) = 0.1, 0.3, 
0.5, 0.7, or 0.9 (from right to left). Our calculations assume a 
constant ionization efficiency in all halos in which atomic cooling 
occurs, set so that reionization ends at z = 6.5. The assumed 
resolution (i.e., the diameter of the region) is 1 or 5 arcminutcs 
as indicated. Note that the range where <5 L < —1 is physically 
sensible since it corresponds to a <5 NL > — 1. 

corresponds in this case to $ max = 0.74), and has most 
of the probability near the cutoff. The reason for this be- 
havior is shown explicitly in Figure [2] Except very early 
in reionization, a small (or even somewhat negative) value 
of overdensity suffices in order to significantly ionize a re- 
gion. In particular, if we increase the overdensity of a region 
in the model, eventually we reach a large-enough overden- 
sity at which the region fully reionizes itself and 9 drops 
to zero. Thus it is not possible to have arbitrarily large val- 
ues of "J. In practice, ^ reaches a maximum at 8 L (z) values 
near zero. Two factors then ensure that most of the prob- 
ability of ^ values is located around this maximum value. 
First, the function ^ versus 8 L (z) is flat near its maximum 
(particularly in the later stages of reionization), and second, 
the probability distribution of S L (z) is centered at values of 
S L (z) near zero. Note that we assume a Gaussian probability 
distribution for S L (z), although the density is weakly non- 
linear and its distribution will thus be slightly modified. For 
instance, the standard deviation of S L (z) in these examples 
is ~ 0.3 for a 1' resolution, and half that for a 5' resolution. 
Since reionization is driven by the distribution of halos, and 
the halo number density is strongly coupled to the mean 
density in each region, we expect the functional form of 
versus S L (z) to be fairly robust. This means that the shape 
of the PDF will also be fairly robust even if the probability 
distribution of density becomes slightly non-Gaussian. 

Returning to Figure [1] we see that in the plotted case, 
the contribution of the one-neutral-one-ionized case to the 
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Figure 3. The CPD of the 21-cm brightness temperature differ- 
ence between two spherical regions. We consider redshift z = 8.26, 
when the mean global ionized fraction is (x l ) = 0.5. We compare 
the full CPD (solid curve), which is the sum of contributions from 
both regions being neutral (short-dashed curve) and from one of 
them being neutral (long-dashed curve), to the CPD of a Gaus- 
sian with the same variance as the total CPD (dotted curve) . The 
comoving separation between the two centers is 10 Mpc, and the 
assumed resolution is 1' or 5', as indicated. In each case, we show 
the purely linear calculation (top panels) or the full calculation 
with a non-linear correction of the density (bottom panels). Note 
that in each panel, the CPD has an additional step function (not 
shown) that brings its value to unity at ATj, = 0. In our cal- 
culations we assume a constant ionization efficiency in all halos 
in which atomic cooling occurs (i.e., in all halos above a circular 
velocity V c = 16.5 km/s), set so that reionization ends at z = 6.5. 



difference PDF is similar (though not identical) in shape to 
the one-point PDF, and in particular is centered near the 
maximum value of AT&. This results from the fact that in 
this case, the value of AT], is simply equal to Tb for the 
region that is not fully ionized. Note that in general the 
maximum value of ATt for the difference PDF is equal to 
the maximum value of Tb for the one-point PDF, since *t is 
always non-negative. The both-neutral contribution to the 
difference PDF is quite different, since $ in each of the two 
regions tends toward $ mal , so the most likely difference be- 
tween the two ^ values is zero. This tendency is further 
strengthened by the correlation between the densities in the 
two regions. The outcome of all of this is a difference PDF 
that has two peaks, with a valley at intermediate values of 
ATb- As a result, the 21-cm difference CPD first declines 
at small values of ATb (where it is dominated by the both- 
neutral case), then flattens at larger values, and finally cuts 
off sharply at the maximum value of AT;,. 

The bottom panels of Figure [3] again show the full CPD 
(with the non- linear correction of the density), but also com- 
pare it to the cumulative of a Gaussian (i.e., an error func- 
tion) with the same variance. The CPD shape described 
above is clearly very different from the error function. As 
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illustrated by the top panels, the non-linear correction that 
we have applied to the density is important in order to en- 
sure that is always positive and that there is a sharp cutoff 
at a maximum value of ATt, (while a calculation with linear 
densities leads to the unphysical result of having negative 
values of $ when the density fluctuation is more negative 
than -1). Other than this cutoff, the non-linear correction 
modifies the shape of the CPD only slightly, with the correc- 
tion having a smaller effect in the case where the resolution 
angle is larger (and where fluctuations on the corresponding 
scale are more linear). Thus, while a non- linear correction 
of the density is required to ensure a physical result (with 
a non- negative density), we do not expect our results to de- 
pend strongly on the precise form of non-linear correction 
that we have used. 

Figure [4] shows the time evolution of the CPD dur- 
ing reionization, considered at two different values of the 
separation d. Throughout the parameter range considered, 
the CPD clearly has the same characteristic shape as noted 
above (although for (x l ) = 0.1 in the top panels, the flat 
portions occur at lower values of the CPD than are included 
in the plot). In the bottom- left panel, we show an example 
of two cases ((a; 1 ) = 0.3 and 0.5) which have nearly identi- 
cal variances (i.e., the corresponding Gaussian CPD curves 
are nearly indistinguishable), and yet the actual 21-cm dif- 
ference CPDs differ substantially in these two cases. The 
Figure illustrates how the CPD evolves during reionization, 
declining with time at low ATi, values (since the probability 
associated with both regions being fully ionized increases), 
and cutting off at a lower value of ATi, (since the overden- 
sity needed for full reionization of a given region declines as 
reionization progresses globally - see Figure [2}. 

The information on spatial correlations contained 
within the CPD is illustrated more clearly in Figure [5] In 
most of the cases shown, the d = 30 Mpc and d = 100 
Mpc curves are nearly indistinguishable, since even at a 30 
Mpc separation the two regions are nearly independent. The 
CPD drops rapidly as the separation is decreased, with the 
probability becoming concentrated near ATi, = once the 
two regions become highly correlated. The decline with sep- 
aration, which occurs at d = 1-10 Mpc early in reionization 
((x l ) = 0.25) but over a broader range of d = 1-30 Mpc 
later on ((x l ) = 0.75), indicates the relative importance of 
bubble and density correlations on various scales. 

A full measurement of the CPD would yield a two- 
dimensional function of d and AT& at each redshift. This 
full function is illustrated with a contour plot in Figure [6] 
Regions that contain much of the probability - i.e., where 
the CPD changes rapidly and there are large spaces between 
consecutive contours - indicate both the characteristic scale 
d of correlations and a corresponding characteristic value 
of ATi, (which is related to the correlated distribution of 
densities in the two separated regions). While the full 21- 
cm difference PDF would be a great tool to study theoret- 
ically and to observe, in the next subsection we show that 
even if the gross features of the PDF were measured, they 
already would reveal important information that is not di- 
rectly available from measurements of the correlation func- 
tion alone. 
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Figure 4. The CPD of the 21-cm brightness temperature differ- 
ence between two spherical regions. We consider redshifts when 
the mean global ionized fraction is {x % ) = 0.1 (dashed curve) or 
(x l ) = 0.3, 0.5, 0.7, or 0.9 (solid curves, from right to left where 
the CPD is small). The comoving separation between the two 
centers is d = 10 Mpc or 30 Mpc, and the assumed resolution is 
1' or 5', as indicated. In the bottom-left panel, we also compare 
the CPD curves to CPDs of Gaussians with the same variance 
(dotted curves, in order from right to left: (x l ) = 0.5, 0.3, 0.7, 
0.1, and 0.9). Note that the CPD has an additional step function 
that brings its value to unity at ATj = 0. In our calculations 
we assume a constant ionization efficiency in all halos in which 
atomic cooling occurs, set so that reionization ends at z = 6.5. 



3.2 Using the difference PDF to separately 
measure ionization statistics 

Next, we consider robust ways to extract information from 
the one-point PDF and the two-point, difference PDF. While 
further study may show that the detailed shape of the PDF 
can be used as a sensitive probe of the underlying astro- 
physics (such as properties of the population of ionizing 
sources), the analytical model we consider in this paper is 
approximate and neglects some non-linear corrections and 
other physical effects, so the precise shape would change 
somewhat in more complete models or numerical simula- 
tions. However, we expect our model to correctly capture 
the gross features of the PDF, which likely constitute more 
robust predictions. Also, on the observational side, while the 
full PDF may be difficult to measure with low signal-to-noise 
data, we expect it to be easier to extract just these gross fea- 
tures. We leave for future study the question of how these 
features can be extracted in a realistic scenario with noise 
(assuming that the PDF of the noise can first be measured 
accurately). Our goal here is to study what information can 
be extracted just from the gross features of the PDF. 

We first consider the one-point PDF. Two gross quanti- 
ties can be simply extracted from it: the probability P(x 1 — 
1) that a region of the resolution size is fully ionized, and the 
value of $ ma x- The first quantity can be extracted by mea- 
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Figure 5. The CPD of the 21-cm brightness temperature differ- 
ence between two spherical regions. We consider redshifts when 
the mean global ionized fraction is (x l ) = 0.25, 0.5, or 0.75, as in- 
dicated. We assume a resolution of 1' (solid curves) or 5' (dashed 
curves), at a comoving separation d = 100, 30, 10, 3, and 1 Mpc 
(from right to left in each set of curves) .Note that the CPD has an 
additional step function that brings its value to unity at AT), = 0. 
In our calculations we assume a constant ionization efficiency in 
all halos in which atomic cooling occurs, set so that reionization 
ends at z = 6.5. 



suring the size of the step function of the CPD, or (equiva- 
lently) by subtracting from unity the value of the CPD just 
above T(, = 0; alternatively, this value equals the integrated 
area under the PDF curve, not including the <5o-function at 
zero. Since most of the probability lies near the maximum 
value of Tt, (see Figure [JJ, measuring it depends mostly on 
pixels with the highest signal-to-noise ratios. Measuring the 
value of ^max also depends on the same pixels; this value 
can be derived from the maximum value of Ti, using eq. 0. 
Note that while a fully realistic PDF may not feature such 
a total sudden drop at a maximum value of Ti,, there should 
nonetheless be a definite, sharp drop due to the strong de- 
pendence of ionization on density through variations in the 
halo number density. Note also that in simulations, while a 
region cannot be truly fully ionized (because of the presence 
of very high-density gas), the interpretation of P(x l = 1) in 
this case is where the bubbles within the region have fully 
overlapped and all the gas is highly ionized except for some 
gas at S ^> 1 (which at high redshift generally makes up only 
a s mall volume and mas s fraction). Indeed, in simulations 
by iMellema et all (|2006l ) the PDF as a function of Ti, has 
a fairly fiat portion (though not always increasing) followed 
by a rather sharp cutoff. We note that the shape of such 
statistics as measured in simulations has not yet been phys- 
ically justified or been subjected to numerical convergence 
tests. But it is generically expected that the two quantities 
P(x l — 1) and 'I'max should clearly feature in the PDF. If 
models or simulations can reliably establish at least the ap- 




AT„ [mK] 

Figure 6. Contour plot of the CPD of the 21-cm brightness tem- 
perature difference between two spherical regions. Same data as 
represented in the previous Figure, but in a contour plot which 
emphasizes the fact that a two-dimensional dataset, Cat 6 as a 
function of AT;, and d, is available to measure at each redshift. 
We again consider (x % ) = 0.25, 0.5, or 0.75, as indicated. Contours 
go from a value C\x b = 1 (short-dashed curve) to Cat,, = 0-01 
(long-dashed curve), passing through Catv =0.1 (dotted curve). 
Equal logarithmic spacing of the values of the CPD is used within 
each decade, with finer spacing at the highest decade of CPD val- 
ues (13 contours total in each panel). We show the case of a 1' 
resolution, and again assume a constant ionization efficiency in 
all halos in which atomic cooling occurs, set so that reionization 
ends at z = 6.5. 



proximate form of the PDF, then this would make it easier 
to measure these two quantities. 

Now consider the difference PDF (see Figure [1} . One 
gross feature is obviously a cutoff that can also be used 
to measure the same value of $ max . The difference PDF 
can also be used to extract three other inter-related quanti- 
ties, each as a function of the separation d: the probability 
P(x\ = 1,352 — 1) of joint full ionization; the probability 
2P(x\ = 1,^2 < 1) of full ionization of only one of the two 
regions (with the factor of 2 accounting for the symmetry 
of selecting which region is the ionized one); and the proba- 
bility P(x\ <l,x\< 1) that neither region is fully ionized. 
The first probability can be extracted from the size of the 
step function of the CPD at ATj, — 0, while the other two 
can be extracted from the areas under the two peaks that 
are fairly well separated in the PDF during reionization (see 
Figure [TJ , or more accurately by modeling the two separate 
contributions to the PDF at ATi, > 0. In reality, the three 
probabilities need not be measured separately, since they 
can easily be shown to be closely related to each other; in 
fact, P(x\ = 1, x\ = 1) together with the one-point quantity 
P(x\ = 1) can be used to express the other two two-point 
probabilities: 

P(x\ <l,xl<l) = 1+P(x\ =l,xl = l)-2P(x l = 1) ,(19) 
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and 

2P{x\ = 1,4 < 1) = 2 [P{x l = 1) - P{x\ = 1,4 = 1)] -(20) 

These relations should make it much easier to extract this 
information from even approximate measurements of the dif- 
ference PDF. 

Thus, from the one-point and two-point PDF, we can 
measure two independent probabilities that depend directly 
only on ionization statistics, not mixed in with the value of 
the density. These are P(x 1 — 1) (a single quantity at each 
redshift) and P{x\ = 1, x\ = 1) (a function of d at each red- 
shift). Even a rough measurement of the PDFs may yield 
a reasonable estimate of these gross quantities. Note also 
that in the limit of infinitely good resolution (i.e., a very 
small resolution scale), P(x 1 = 1) becomes the cosmic mean 
ionized fraction, and P(x\ = l,x\ = 1) becomes the ioniza- 
tion correlation function (after subtracting the square of the 
mean ionized fraction). In addition, the value of $ max yields 
an interesting piece of information on the dependence of ion- 
ization on the density. All of these quantities are separate 
from the correlation function, which yields just one func- 
tion of d that is a complicated convolution of fluctuations in 
density and ionization. 

Figure [7] shows predictions of our model for the two 
quantities P(x l = 1) and *l/ m ax that are available from the 
1-pt PDF. P(x l — 1) can be used as a rough estimate for 
the cosmic mean (x l ), although this works better with high 
resolution (1') and late in reionization. The robustness of 
theoretical predictions of the relation between P(x z = 1) 
and (x z ) can be investigated with further models and sim- 
ulations. Even for a fixed end-of-reionization redshift and 
when expressed as functions of (x z ), both P(x l = 1) and 
*Pmax depend significantly on the characteristic halo mass of 
ionizing sources. 

Figure [8] shows the three inter-related probabilities ob- 
tainable from the difference PDF, compared to the expec- 
tation value ($1^2)- The standard 21-cm 2-pt correlation 
function does not actually yield {^1^2) but rather subtracts 
off (*) 2 , which is the value of (*i* 2 ) at d —> 00. Thus, in 
the figures the information available from the 2-pt correla- 
tion function is not the absolute plotted values of {9 1^2), 
but just the values relative to the large-scale asymptotic 
value. Therefore, the ability to measure the 2-pt correlation 
function and extract useful information from it depends on 
the total change in (^1^2) from small to large scales. This 
total change is smaller than the overall change with scale 
of most of the curves that show the ionization probabilities. 
The probability P(x\ = l,x 2 = 1) goes from P(x\ = 1) at 
d -> to [P{x\ = l)] 2 at d -> 00, while P(x\ < 1,4 < 1) 
goes from 1 - P(x{ = 1) to [1 - P{x\ = l)] 2 . The largest 
change is seen in the probability that one region is ionized 
and the other is not; this varies from zero at small d to 
2{P(x\ = 1) - [P{x\ = l)] 2 } at large d. 

The characteristic scales of the bubble correlations can 
be read off as the scales where the various quantities in 
Figure [8] change most rapidly as a function of d. In or- 
der to focus on this important feature, we plot the deriva- 
tives with respect to log(d) in Figure [9] just for {^1^2) 
and for P(x\ = l,x 2 = 1) (since the other two probabil- 
ities can be derived from this quantity). In general, both 
P(x\ = l,x\ = 1) and (^1^2) show roughly the same dom- 
inant scales, but P(x\ = 1, x\ = 1) shows a greater variation 
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Figure 7. Separate observables available from the 21-cm one- 
point PDF. Assuming parameters for which the universe fully 
rcionizes at 2 = 6.5, we consider stars forming (with a constant 
ionization efficiency) in all halos above V c = 16.5 km/s (corre- 
sponds to efficient atomic cooling) or only in halos above V c = 45 
km/s (corresponds to strong feedback in low-mass halos, e.g., due 
to photoheating or supernovae). In each panel, we show P(x l = 1) 
for V c = 16.5 km/s (solid curve) or V c = 45 km/s (dot-dashed 
curve), compared to the cosmic mean {x l ) (straight dotted line). 
We also show ^max for V c = 16.5 km/s (short-dashed curve) or 
V c = 45 km/s (long-dashed curve). 

with scale during the central stages of reionization at which 
the variation is maximized. Figure[lO]shows similar trends in 
the case of reionization by more massive and highly-biased 
halos, except that here the characteristic scale grows to even 
larger values by the end of reionization; early on in reioniza- 
tion, the characteristic scale is approximately the same as in 
the previous case of less-massive halos, but the magnitudes 
of the scale-derivatives are larger in the case of the more- 
massive halos (corresponding to stronger ionization and 21- 
cm fluctuations). Figure [TT1 demonstrates the importance of 
achieving high resolution in the upcoming experiments. In 
the scenario considered here, a cosmic mean ionization frac- 
tion of one half occurs at z = 8.26, at which the resolution of 
V corresponds to a radius of 1.3 comoving Mpc (com Mpc), 
and 5' to 6.7 com Mpc. The typical correlation scale grows 
with time and overtakes the 5' scale only late in reioniza- 
tion ({x" 1 ) ~ 0.7), so until then the peaks of the curves in 
Figure [TT] indicate roughly the resolution scale (rather than 
the desired correlation scale). Observations should therefore 
achieve a resolution of a few arcminutes or better in order to 
make it possible to measure the evolution of the dominant 
correlation scale throughout most of the reionization era. 



4 CONCLUSIONS 

We have presented and studied a new statistic for analyz- 
ing 21-cm fluctuations, namely the PDF of the difference 
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Figure 8. Gross observables available from the 21-cm two-point 
difference PDF, shown versus the separation of the two points. 
Assuming parameters for which the universe fully reionizes at 
z = 6.5, we consider stars forming with a constant ionization ef- 
ficiency in all halos above V c = 16.5 km/s. For the resolution 
we assume a value of 1'. We consider at various stages dur- 
ing reionization {(x l ) = 0.1, 0.3, 0.5, 0.7, or 0.9), the standard 
quantity (\Pi\E r 2) from the 21-cm two-point correlation function 
(dotted curves, from top to bottom) as well as the three inter- 
related gross quantities from the difference PDF (which together 
add up to unity): P(x\ = l,x l 2 = 1) (solid curves, from bot- 
tom to top), 2P(x\ = X,x\ < -0 (long-dashed curves; at d=100 
Mpc, (x») = 0.1, 0.9, 0.3, 0.5, and 0.7, from bottom to top), and 
P(x\ < 1,^2 < 1) (short-dashed curves, from top to bottom). 

of the 21-cm brightness temperatures of two regions, as a 
function of the separation between their centers. This two- 
dimensional statistics generalizes to one-dimensional func- 
tions, the one-point PDF and the two-point correlation func- 
tion; the latter is simply related to the variance of the dif- 
ference PDF (Eq.QSJ. 

We have predicted the difference PDF based on the 
correlated two- point distribution of density and ionization 
l|Barkanai r2007). generalized to the case of a finite observed 
resolution, and including a non- linear correction for the den- 
sity (Eq. [5} . The model predicts a characteristic shape dur- 
ing reionization for the PDF (or its cumulative form, the 
CPD) that can be understood from the various contribu- 
tions to it depicted in Figure [T] The PDF contains informa- 
tion on the distribution of the density within neutral regions 
and on the dominant spatial scales of bubble and density 
correlations (Figures U and [SJ. The full PDF is a function of 
separation and of AT;, at each redshift (Figure [6}. 

While the usual correlation function is determined by a 
complicated mixture of contributions from density and ion- 
ization fluctuations, we have shown that the difference PDF 
(together with the one-point PDF) holds the key to sepa- 
rately measuring statistics of the ionization distribution. In 
particular, even an approximately measurement of the PDFs 
can generically be used to measure the ionization probability 
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Figure 9. Derivatives of the 21-cm two-point quantities. We 
show at various stages during reionization \d / d\og(d)\ of P(x\ = 
l,x 2 = 1) (solid curves) and of (^l*^) (dotted curves). For 
each quantity, we consider various stages during reionization 
({x 1 ) = 0.1, 0.3, 0.5, 0.7, and 0.9, in order from left to right 
at the peak of the curve). 

of a resolution-sized region, and the joint ionization proba- 
bility of two such regions as a function of their separation. 
Within our model, the joint ionization probability shows the 
same characteristic correlation scale as the two-point 21-cm 
correlation function but has more power (Figures l9l and [T0|) : 
this is because the contributions of density and ionization to 
21-cm fluctuations are anticorrelated (a higher density im- 
plies less neutral gas), which reduces the 21-cm fluctuations 
relative to ionization fluctuations. 

If quasars contribute significantly to reionization by 
producing large bubbles even early in reionization, or if 
quasars or supernovae emit significant X-ray photons (which 
have a long mean free path), then the density and ioniza- 
tion fluctuations will not be as simply related as we have 
assumed. In this case, it will be much more important to 
measure the ionization statistics separately since it will be 
difficult to extract this information from the 21-cm statis- 
tics. Also in this case the relation between the one-point and 
difference PDFs and the 21-cm correlation function should 
be significantly different from the case we have considered. 
The full quantitative details of such a scenario goes beyond 
the scope of this paper and merits a separate study. 
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Figure 10. Same as Figure [9] except that we assume stars form 
with a constant ionization efficiency in all halos above V c = 45 
km/s. 
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Figure 11. Same as Figure[9] but for the case of a 5' resolution. In 
this case, the (x l ) = 0.1 case cannot be seen in P(x\ = 1, x l 2 = 1) 
(which is essentially zero in this case, at all d). Here, the order 
of increasing peak heights (not positions) is (x % ) = 0.1, 0.3, 0.9, 
0.5 and 0.7, for (*i* 2 >, and (a:*) = 0.3, 0.5, 0.9 and 0.7 for 
P(x\ = 1,4 = 1). 
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